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Abstract. We discuss the laser theory for a single-mode laser with nonlinear gain. 
We focus in particular on a micromaser which is pumped with a dilute beam of 
excited atoms crossing the laser cavity. In the weak-coupling regime, an expansion 
' in the coupling strength is developed that preserves the Lindblad form of the master 

equation, securing the positivity of the density matrix. This expansion breaks rapidly 
CN \ down above threshold. This can be improved with an alternative approach, not 

restricted to weak coupling: the Lindblad operators are expanded in orthogonal 
. polynomials adapted to the probability distribution for the atom-laser interaction time. 

Results for the photon statistics and the laser linewidth illustrate the theory, 
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1. Introduction 

The quantum theory of a laser is a textbook example of a nonlinear problem that 
requires techniques from open quantum systems. The key issue is the nonlinearity 
in the gain of the laser medium, due to saturation, that leads to coupled nonlinear 
equations already at the semiclassical level. The quantum theory makes things worse 
by its use of non-commuting operators. 

Recall that in the so-called semiclassical theory, the following equation of motion 
for the intensity / the laser mode can be derived flU [2D : 

d/ r AI 

where k is the loss rate, A is the linear gain, and (3 describes gain saturation for the 
laser medium. A quantum upgrade of this theory replaces the intensity by the photon 
number a) a where the annihilation operator a describes the field amplitude of the 
laser mode. Mode loss is easy to handle by coupling the laser mode linearly to a 
mode continuum 'outside' the laser cavity [3]. This leads to a master equation for 
the density matrix in so-called Lindblad form [see Eq.©] with a Lindblad operator 
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L = y/R a. Linear gain can be handled in the same way, but gain saturation is more 
tricky. A heuristic conjecture is a Lindblad operator L = y/Aa)(l + /3a^a) -1 / 2 . The 
operator ordering can only be ascertained a posteriori, and it is difficult to choose 
among the replacements / i— > a^a, aa^, or \{a)a + aa)}. 

The pumping of the laser can be modelled in different ways, for example by 
placing excited two-level atoms into the laser cavity. Nonlinear gain emerges from 
a treatment beyond second order in the atom-field coupling. In Ref.[4], a coupling 
to fourth order is considered, and in addition, an average over the atomic lifetime 
in the cavity is performed. In Refs.Q2l [5D, the pumping model is based on a dilute 
stream of excited two-level atoms that cross the laser cavity one by one and interact 
with the laser mode during some randomly distributed interaction time. This model 
can be largely handled exactly (61 , even in the presence of incoherent effects like 
cavity damping, imperfect atom preparation, and frequency-shifting collisions. The 
setup has become known as the 'micromaser' because of its experimental realization 
with a high-quality cavity (71 [8l [9j] . One line of research has focused on the so-called 
'strong coupling regime' that permits the laser mode to be driven into non-classical 

states mum. 

We focus here on the 'weak coupling' regime. On the level of the master equation 
for the laser mode, this regime corresponds to a small product of coupling constant 
and elementary interaction time r so that one can expand in this parameter. For the 
description of a realistic experiment, one has to average the master equation with 
respect to a distribution in r (Sec. [2]) . It turns out, however, that the resulting master 
equation is not of the well-known Lindblad form, although it preserves the trace of 
the density matrix. This leads to conflicts with the positivity of the density operator, 
as is known since the original derivation of the master equation by Lindblad and by 
Gorini et al. (121 [T3Tl . In this paper, we give a discussion of this problem and suggest 
a solution. On the way, we review the derivation of the Lindblad master equation 
starting from the Kraus-Stinespring representation of the finite-time evolution of the 
density matrix (Sec. [3]). The mathematical treatment is at the border of validity of 
the formal Lindblad theory since one has to deal with an infinite-dimensional Hilbert 
space and continuous sets of Kraus and Lindblad operators. We construct two modified 
expansions that result both in a Lindblad master equation (Sees. 01 [5]). One is a 
direct amendment of the weak coupling approximation, the other one is able to enter 
the regime of a strong coupling (on average). The latter gives at least qualitative 
agreement with the results of the exact master equation. 

2. The micromaser model 

Consider a two-level atom with states \g), |e) that is prepared at time t in its excited 
state |e) = (1,0) T (density matrix pa = |e)(e|) and that interacts with a single mode 
(density matrix p) during a time r. One adopts a Jaynes-Cummings-Paul Hamiltonian 
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for the atom-field coupling 



)■ 



a = \g)(e 





1 



) 



(2) 



(this applies at resonance in a suitable interaction picture). Assume that the initial 
density operator of the atom + field-system factorizes into P(t) = pit) ® p A , compute 
P(t + t) by solving the Schrodinger equation and get the following reduced field 
density matrix fl2l [5]] 



where sinc(x) = sin(x)/x, and (p 2 = aa) is one plus the photon number operator. The 
operator-valued functions cos and sine are defined by their series expansion. Only 
even powers of the argument occur, hence we actually never face the square root (p of 
the operator aa 1 ". In the following, we abbreviate the mapping defined by Eq.© by 
M T p(t) (this is sometimes called a superoperator) . 

The operation O) describes an elementary 'pumping event' of the laser. To 
provide a more realistic description, one introduces the following additional averages: 
excited atoms appear in the laser cavity at a rate r such that rr < 1. The interaction 
time t is itself distributed according to the probability measure dp(r) with mean value 
f . On a coarse-grained time scale At 3> f , this leads to the difference equation [j2l ED 



To simplify the superoperator appearing on the right hand side, Refs. []2l |5B suggest an 
expansion in powers of gr(p up to the fourth order. Using an exponential distribution 
for dp(r), this leads to the approximate master equation 



where we followed the common practice of interpreting this as a differential equation. 
We use {•, •} to denote the anticommutator. The linear gain is A = 2r(gf) 2 , and 
B = (gf) 2 A is a measure of gain saturation. Losses from the laser mode can be 
included in the usual way by adding a term of the same structure as the first line of 
Eq.([5]), but exchanging a and and replacing A by the cavity decay rate fl2l [5D . The 
same master equation is also found, using a different pumping model, in Ref. fl4D . 

It is easy to check that Eq.© preserves the trace of p, using cyclic permutations. 
Nevertheless, it is not of the general form derived by Lindblad for master equations 
that preserve the complete positivity of density matrices Ifl2[ [T3l [14ll . We shall 
show below fEq. (f36l) 1 that Eq.© indeed leads to a density matrix with negative 
probabilities. Recall that the Lindblad form is given by 



p(t + r) = cos(gr<p)p(t) cos(gnp) + (gr) 2 a 



smc(gT0)p(t) smc(gT0)a (3) 




(4) 




(5) 




(6) 
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with a countable set of operators L x . One may think of an ansatz polynomial in the a 
and a) for the L x , but it is difficult to see how to generate the mixed third order-first 
order terms a^aa^pa in Eq.© without generating also contributions like cJaa) paa)a. 
Note that this 'missing term' cannot disappear by cancellations: if we allow the L x 
operators to contain at maximum three factors of a or then the highest order term 
generated by the 'sandwich' structure L x ph\ is proportional to the squared coefficient 
of the highest order term of L x , and these terms cannot cancel out. 

Of course, one can accept to work with this kind of 'post-Lindblad' master 
equations (as they appear frequently in the papers of Golubev and co-workers, see 
e.g. lfT5lO . We follow here another route and raise the question: what assumptions 
behind the standard Lindblad master equation do not apply here, or are we missing 
something? To formulate an answer, we go back to a derivation of the Lindblad 
form that starts from another general formulation for mappings between density 
matrices, the so-called Kraus or Stinespring representation of completely positive 
operators [14]. We show that a set of Lindblad operators {L x } can indeed be 
constructed so that by adding a few additional terms to the master equation ©, it 
can be brought into the Lindblad form. 

3. Lindblad from Kraus-Stinespring 

The time evolution of a density matrix can be expected to yield a density matrix again. 
This intuitively obvious requirement is violated by some models 111611 or for some 
initial states |fT7lK but it can also be taken as a starting point for an 'axiomatic' theory 
of dissipative quantum dynamics. Following this latter approach, one derives that the 
evolution over a finite time At must be of the form (Stinespring theorem, Ref. HT4l ) 



where the operators Q x depend on At and satisfy the 'completeness relation' 
Xa^a^a = 1 to ensure trace conservation. This form (called Kraus representation 
11141 [1811 ) can be easily secured for the micromaser master equation ([4]). We resolve 
the discrete difference quotient and get (A = 0, 1, 2) 



where the completeness relation is satisfied because of the trigonometric identity 
sin 2 + cos 2 = 1 that is carried over to operator-valued arguments. 

We note that the Kraus representation retains its form, at least formally, when we 
average the operators f2 A with respect to a distribution in the parameter r. This is 
easily seen by interpreting the integral over rasa Riemann sum: for each A, the term 



p(t+At) = Y,n xP m\ 



(7) 



A 



tt = (1 -r At) 1/2 1 

fti = (rAt) 1/2 cos(5fT<£) 

^2 = (rAty^gra^ sinc(grip) 



(8) 
(9) 
(10) 
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Q x (r)pQ x (r) is replaced by the sum 

E^V^l- with % = ^A(r i )V / dK^) (ID 
j 

The completeness relation is also still satisfied: for each tj, f4( T i)^i( T j) and 
Ql(Tj)Q 2 ( T j) stn l a dd up to the unit operator, and the sum over the prefactors goes 
over into the normalization integral of the probability measure dp(r). Note that we 
interchange here the summations over A and j. 

The derivation of the Lindblad master equation (see, for example, Ref.[18] and 
the Appendix) now provides a construction of the operators L\ appearing in Eq.©. 
Extract a traceless operator V\ by writing 

Q x = lj x 1 + V x (12) 

and take the limit 

L x = lim -^=. (13) 



The operator fi is already proportional to the unit matrix, hence Vo = 0. It is also 
obvious that tt 2 is traceless, hence 

L 2 = S = y/r gra * sxnc(gr<p) (14) 

The operator Q 1 has a singular trace. In the number state basis: 

oo oo 

J2H^i(r)\n) = (rAt) 1/2 ^ cos(prVn + l) (15) 

n=0 n=0 

The subset of square numbers gives a divergent result whenever gr = modulo 2-n, 
hence a 'comb' of 5-functions is expected. This is of course a tricky result in view of 
the expansion in gr that is operated in the way from Eq.© to Eq.©. We therefore 
introduce a factor q n with < q < 1 into the sums (fT5l) . Comparing the traces of both 
sides in Eq. (fT2~l) for A = 1 (using that V\ is traceless), we get 

oo 

wi = (rAt) 1/2 w(gr) = (rAt) 1/2 ]T (1 - q)q n cos^rVnTT). (16) 

ra=0 

Finally we find 

L x = C = y/r [cos(gnp) - zu^gr)!} (17) 

where zui(gr) is defined in (fT6l) . 

Observe that at this stage, we do get a master equation in Lindblad form. But the 
Lindblad operators still contain the interaction time gr to all orders. 



4. Weak-coupling Lindblad form 



We now investigate how the expansion in powers of gr and the averaging with respect 
to dp(r) can be organized so that the Lindblad form is preserved. 
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4.1. Consistency of the expansion 

We start with two general remarks. Consider a polynomial approximation of order N 
to the cos and sine functions in Eqs. (fT4l[T7l) . The operator C^C is then of order 2N in 
gr and the operator S^S of order 2N + 2 [cf. Eqs. ([T4lfT7l )l. The maximum number of 
factors a or a* in the master equation is given by 2N + 2. A scheme consistent with the 
Lindblad form thus seems possible only if the master equation is expanded at least to 
the order 2N + 2 in gr. 

The case 2N + 2 = 2, involves the second order (gr) 2 only, hence a Lindblad-like 
form with rate coefficient 2r(gf) 2 = A. This reproduces the first line of the master 
equation d5]). 

The next case is 2N + 2 = 6 because both cos and sine are even in r. Then one 
should have six factors a or a) . We see that this is not the case in Eq.© where only 
up to four factors appear. Therefore, the expansion in gr(p has not been pushed to 
a high enough order (sixth) to be compatible with a Lindblad form, at least for the 
terms originating from the Lindblad operator S. 

A second remark: consider the expansion in powers of gr of the operators C, S: 

oo 

C = ^Y,(9T) n C n (18) 

n=0 

oo 

S =Vt^J2(9T) n S n (19) 

n=0 

where the coefficients C n ,S n involve powers of the operator <p = (aa 1 ") 1 / 2 . The 
integration over r now leads to 'cross terms' like 

Jdp{r){gr) n+m C n pC m . (20) 

The resulting master equation is not in diagonal form because once the integration 
over t performed, these cross terms cannot be written as a product of a function of 
n times a function of m. Progress can be made by using an expansion in orthogonal 
polynomials, as we discuss now. 



4.2. Polynomial expansion 

In the expansions (fTSl [T9l ) of the Liouville operators, let us re-write the powers as 

n 

r n = J2 a nkf n f k (r/f). (21) 

k=0 

The factor f n is chosen for dimensional convenience. In the example discussed below, 
f is identified with the mean value of the probability measure dp(r). The polynomials 
fk(r/f) are of order k and are orthogonal with respect to the following scalar product 

dp(r)f k (r/f)f l (r/f) = 6 kl . (22) 

This is a scalar product since dp(r) is a positive measure. Such polynomials exist and 
are real. An explicit example is worked out below for an exponential distribution. The 
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coefficients a nk in Eq. (f2~T1) can be found by projecting x 11 onto f k (x) which boils down 
to an integral similar to (l22l) . Note that a nk = for k > n because x n can be written as 
a finite linear combination of the f(x) (0 < I < n), using the Gram-Schmidt procedure 
for orthogonalization. We then get for the first cross term in (l20l) : 

min(n,m) 

dp(r)r n+m = f n+m J2 a nta mk . (23) 

k=0 

The average of the master equation involving the Lindblad operator C, say, then 
assumes the following diagonal form 

J dp(r) (Cptf - \ {tfC, p}) 

oo oo 

= rJ2 E *n k a mk (gf) n+m (C nP Cl - 1 {Clc nl p}) (24) 

fe=0 n,m=k 

which can be written in a Lindblad form involving the (countable set of) operators 

oo 

C k = ^T. a nk{gr) n C n . (25) 

n= k 

If the expansion (fT9l) is truncated at order N, then C k involves also only terms up to 
order n = N and the set of Lindblad operators is finite as well. Let us consider N = 2 
and take into account that C is even in r. Then 

Co = aooCo + a 20 (gf) 2 C 2 (26) 
Ci = a 2l {gffC 2 (27) 
C 2 = a 22 (gf) 2 C 2 (28) 

The Lindblad operators C\ 2 can be combined into a single one since they are 
proportional to the same operator C 2 . A similar procedure can be applied to S, the 
only difference being that only odd coefficients Si, S3, ... are nonzero. 

4.3. Example: Laguerre polynomials 

The Laguerre polynomials L n (x) implement orthogonality with respect to a scalar 
product weighted with an exponential |fl9ll 

00 

J dxe~ x L n (x)L m (x) oc S nm (29) 


which corresponds to the probability distribution dp(r) = (dr/f )e~ T//r considered by 
Mandel and Wolf [4] and by Orszag [2J. We identify x = r/f as the natural variable 
for the polynomials we require. The first few Laguerre polynomials read, normalized 
as in Eq.d22]) 

f (x) = l fi(x) = l-x 

f 2 ( x = ( x 2 - Ax + 2 A 
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A straightforward calculation gives the following Lindblad operators for the master 
equation {[5]): 

Co = M9r) 2 (jYZ-qj - aflt ) ' ^ = -2^0, C 2 = C (31) 
S = ^gfa ] ( 1 - (gffaa)) (32) 
§!= - yfigfa) ( 1 - 3(gf ) 2 aa t ) (33) 
S 2 = vToK#r) 3 aW (34) 

The operators Co,i,2 are proportional to each other and can be combined into a single 
one (replace ^fr by y/Qr in (%). An analoguous simplification has been already made 
in writing Eq. (f34l . Working out the details, we see that the part of C that involves 
1/(1 — q) actually does not contribute to the master equation. (This is generally true 
if we have a hermitean Lindblad operator and add a term proportional to the unit 
operator with a real coefficient.) 

We can now identify the 'missing pieces' in the master equation fl5]) . Collecting 
the third-order terms arising from £0,1,2 gives 

-f = 20 r(gff (0+00+ p acta - §{(aa f ) 3 , p}) . (35) 

°t 6th 

These terms are, of course, of sixth order in (gf) 6 and, not really surprisingly, 
themselves in Lindblad form. All other terms are of lower order in gf and combine to 
reproduce Eq.©. 

4.4. Numerical results 

To illustrate the accuracy of the expansion performed here, we have worked out 
the equilibrium photon statistics, i.e., the diagonal elements p n = (n\p eq \n) of the 
stationary solution to the master equation. Two examples are shown in Fig. [H for 
the same value of the pumping parameter A/k and different values of the coupling 
strength gf. The photon statistics is fairly well approximated at weak coupling 
(gf = 0.03), as expected. At the value gf = 0.15, the average photon number is 
not very large, and significant differences occur. 

To understand these differences, consider the following recurrence relation that 
determines the photon statistics in the case of the weak-coupling approximation (see 
Refs.EUHI!). 

Pn+1 = 2 -^- (l - A{gf f{n + 1) + lOfof )\n + l) 2 ) Pn (36) 

Note that the factor in parentheses is positive for all n and becomes larger than unity 
for n 3> n cnt = ^(gf)~ 2 , leading to a divergence of p n at large photon numbers. To 
enforce convergence, we have cut off the number distribution at n cut . This does not 
change the results if this number is well beyond the peak of p n (weak coupling) . But 
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Figure 1. Equilibrium photon statistics for the micromaser. Thick solid line: exact 
theory Q, as worked out in Ref.[2]. Dashed line: Lindblad theory in the weak 
coupling limit with operators G31H34D . Thin solid line: Lindblad theory in the uniform 
approximation ( 1401142 1 ) , see Section [5j In the right plot, the symbols mark the discrete 
values for the photon number. The average interaction time is fixed to gf — 0.03 (left) 
and gf — 0.15 (right). 



as gf increases, the probabilities p n (n pa n cut ) near the cutoff are still significant, and 
the approximation breaks down. 

The photon statistics also allows us to illustrate the failure of the non-Lindblad 
master equation dH). This theory leads to a recurrence relation identical to Eq. (f36l) . 
except that the last term in the parenthesis is missing. This leads to negative 
probabilibites p n for n > \(gf)~ 2 . This obviously unphysical result is a clear 
manifestation of a non-positive density operator, while the preservation of positivity 
is a key assumption in the derivation of the Lindblad form. In Ref.[4], this problem is 
circumvented replacing Eq. (l36l) (without the last term) by 



which does not violate positivity. Incidentally, the description then becomes equivalent 
to the Sargent-Scully laser theory [fill . 

The average photon number (n) and its normalized variance Q = (An) 2 / (n) 
(essentially the so-called Mandel parameter) are plotted in Figs. [2] and [3l A 
similar trend can be observed, with the weak coupling expansion giving an accurate 
description below and slightly above threshold. The agreement is the better, the 
smaller the coupling parameter gf. Above threshold, the expansion is no longer useful 
because photon numbers with gfy/n + 1 ~ 1 are significantly populated. At threshold, 
the photon number fluctuations are strongly super-Poissonian (the Mandel parameter 
Q > 1). They tend to the Poisson (or coherent state) limit above threshold, but this 
regime is not accessible with the weak coupling expansion. We develop an alternative 
description (leading to the thin solid lines) in the following Section. 

Finally, we plot in Fig|4] the following estimate for the laser linewidth 




(37) 




(38) 
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Figure 2. Average photon number (i.e., laser output intensity) vs. pumping strength. 
Left: weak coupling, gf = 0.03; right: stronger coupling gf = 0.15. 
Thick solid line: exact theory; dashed line: Lindblad theory for weak coupling; thin 
solid line: uniform Lindblad theory. 




0.5 1 1.5 2 0.5 1 1.5 2 2.5 3 3.5 



Figure 3. Normalized variance of the photon number, Q = (An) 2 /(n) (Mandel 
parameter), vs. pumping strength. The curves are labelled as in Fig|2l 

where the derivative with respect to r is evaluated using the master equation and we 
consider t — > oo so that the expectation value is taken with respect to the stationary 
state. The data plotted in the Figure are normalized with respect to k/(ti) which 
is of the order of the Schawlow-Townes linewidth. Values close to one indicate the 
line narrowing typical for a laser above threshold. We see that the weak coupling 
approximation rapidly deviates above threshold. At strong coupling, significant 
deviations from the Schawlow-Townes limit occur in all descriptions. This can be 
traced back to an additional, positive contribution from to the C-operators in the 
master equation. Note that in both the exact and approximated theory, these operators 
are diagonal in the number state basis and hence do not influence the photon statistics. 
We shall report on a more detailed analysis elsewhere. 

5. Uniform expansion 

To conclude, we discuss an alternative expansion for the Lindblad operators. The 
idea is to perform an expansion of the operators C{gr(p) and S(grtp) in Laguerre 
polynomials in r. The average with respect to dp(r) is then easy due to the 
orthogonality relation (|29l) . This gives a different dependence on <p and gf where 
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Figure 4. Normalized linewidth D of the micromaser, vs. the pumping strength. The 
linewidth is normalized to the Schawlow-Townes limit n/(n). The curves are labelled 
as in FiglJl 



actually the operator (p appears to all orders. We shall see that this approximation 
provides a convergent photon statistics even above threshold. 

For simplicity, we retain in the expansion only the lowest order polynomials and 
approximate the operator S by / (r/f) ^.uniform + /i(T/f)Si )Um f orm . We have 

oo 

Sk,nniiorm = J dp(r) S (gnp) f k (r/f), A; = 0,1, (39) 



which results in 



^uniform = y/f gT CL J r, (40) 

1 + {gT) z aai 

5 r - t 1 ~ (^) 2 QQ f r „-n 

^ 1, uniform = ~ \J T gT a 1 -r-. (41) 

[1 + (gTyaai\ z 



Note that both operators contribute at order V A for small photon numbers such 
that (gf) 2 (n + 1) <C 1. The reduction of the weak signal gain occurs through 
the denominators that involve the photon number, similar to the Scully-Lamb laser 
theory [ 1 ] . But observe that gain saturation even happens with the laser mode in the 
vacuum state. This regime has been studied previously to prepare, e.g., non-classical 
states in the micromaser. In the visible frequency band, the regime is accessible for 
microlasers with high-Q cavities [120), [2TH . 
A similar calculation leads to 



a 



0, uniform 



I V v v ^ (42) 



which features the same gain saturation. Observe again that the part proportional to 
the unit operator actually drops out from the master equation. The following orders 
(involving the polynomials fi^(x)) are proportional to ^Jr(gf) 2 = y/A(gf) and can be 
included to systematically improve the approximation. 

The matrix elements of these 'uniform' operators are decreasing as the photon 
number gets large. This provides a recurrence relation for the photon statistics that 
converges for n — > oo. We indeed observe from the numerical results shown in 
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Figs. [TH3] that the divergences of the weak coupling approximation are removed. The 
behaviour of the exact theory for all considered observables is well reproduced. We 
speculate that additional terms in the Laguerre expansion will improve the results for 
the laser linewidth (Fig. [4]) where the agreement is worse than for the moments of the 
photon statistics (Figs. El [3]). 

Let us finally comment on the heuristic choice L = \f~Aa\l +f3a^a)~ 1 ^ 2 mentioned 
in the Introduction. The choice f3 = 4(gf) 2 leads, by construction, to the same 
photon statistics as the exact theory. It does not feature the onset of gain saturation 
already for the vacuum state (as the micromaser theory does), unless one changes the 
order of operators. Another shortcoming is the laser linewidth that is not correctly 
reproduced, as additional contributions arise from the C-type Lindblad operators. 
Hence, for the micromaser at hand, this approximation is not suitable. It can be 
used as an introductory tool for more conventional lasers, with the advantage that one 
automatically gets a master equation that is trace preserving and whose rate equations 
satisfy detailed balance. 

6. Conclusion 

Lasers with a nonlinear gain are typically modelled by coupling a reservoir of excited 
two-level atoms to the laser cavity. A specific realization is the micromaser where 
a dilute jet of excited two-level atoms crosses the cavity. This has been studied 
both experimentally and theoretically for a long time already. We have pointed out 
here that an expansion of the master equation in the weak coupling regime can 
be organized in such a way that the equation retains its Lindblad form explicitly. 
This automatically avoids unphysical predictions involving, for example, negative 
probabilities. The expansion is based on polynomials that are orthogonal with respect 
to the probability distribution of the atomic transit time through the laser mode. An 
alternative scheme that is able to handle the strong coupling regime as well has been 
suggested and leads to a reasonable agreement with the exact theory. Further work 
will address a detailed analysis of the laser linewidth and the strong coupling regime. 

Acknowledgements. We thank G. Morigi for helpful comments on the manuscript, 
and M. Martin and T. Felbinger for hints on the numerical calculation. Constructive 
remarks from the referees are gratefully acknowledged. The Appendix is adapted from 
Chap. 5 of Ref.H22E. 

Appendix A. Derivation of the Lindblad form 

In an axiomatic approach to the time evolution of density matrices, it can be shown 
that over a time At, the density matrix changes according to Eq. d71) lTT4l [T8l1 . The 
Lindblad theorem then states: 



Laser theory in manifest Lindblad form 13 
Suppose that the time evolved density operator has the weak continuity property 



lim 

At-»0 



Ap(t + At) - Ap(t) =0(At) (A.1) 



for all operators A and initial density matrices pit). Then there exists a 
hermitean operator H and a set of traceless operators L\ such that 

d P :rrr _i , Wrt„r Ut rt 



= -i [H, p]+E [L\pL x - |{L A Ll, p}) (A.2) 

This differential equation is called the Lindblad form and the L x are called 
Lindblad operators. 

Proof. Let At > and write p = p(t) for simplicity. We start with the Kraus 
representation (0) for the density matrix pit + At), 

p(t + At)=Y.n xP n{ (A3) 

A 

The operators occurring in Eq. (IA.3D can be split into 

Q x = u x l + V x (A.4) 

where the V\ are uniquely defined by the requirement that their trace be zero. Note 
that io x and V\ depend in general on At. 

In terms of these quantities, the change in the density matrix is computed to be 



p(t + At)-p = (^2 -iJp + E (colV xP + 



pU X Vy 



+ £ VIpVx (A.5) 

A 

where u x is complex conjugate to co x . Using the continuity condition (LAID for all 
operators A and p, we find 

lim£> A | a = l (A.6) 



At^O ' 

A 



lim Y oo* x V x = (A.7) 

A 

hE^^ = (A.8) 



At->0 , 
A 

where the last line applies to any density matrix p. We can thus introduce the 
derivatives 

7 = lim ^ A ^ 2 - 1 (A.9) 

' Ai-^0 At 

Y-iH= lim a (A.10) 
At^O At 

where T and H are both hermitean. 
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Differentiating the condition that the dynamical map preserves the trace of the 
density matrix, we find 

tr [p(t + At) - p] 



= lim 



At->0 At 
tr 



IP + 2I> + Km i Vfc p] (A. 1 1) 

A 

Since this must hold for any density matrix p, we find another derivative 

lim ^ A A A A = -7 - 2Y (A.12) 

Ai^O At ' 

We can thus introduce the Lindblad operators L A by the limiting procedure 

L\ = lim -fL= (A. 13) 

a*-»o VaT 



Using the derivatives defined in Eqs. flA.9[ IA.101 IA.13D , we can divide the difference 
p(t + At)-p(t) in Eq. (lA.5D by At, and take the limit At -> 0. This gives the differential 
equation OA. 2D . ■ 
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